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ABSTRACT 

A  new  iterative  method  is  presented  for  solving  finite  difference 
equations  which  approximate  the  steady  Stokes  equations*  the  method  is  an 
extension  of  successive-over-relaxation  and  has  two  iteration  parameters* 
Perturbation  methods  are  used  to  analyse  the  iteration  matrix.  Sufficient 
conditions  for  the  convergence  of  the  iterative  method  are  obtained  and  it  is 
shown  that  many  reasonable  finite  difference  schemes  for  the  Stokes  equations 
satisfy  these  conditions.  Computational  examples  are  given  to  show  the 
efficiency  of  the  method. 
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SIGNIFICANCE  AND  EXPLANATION 


The  incompressible  Navier-Stokes  equations  describe  the  flow  of  many 
common  fluids.  Thus  effective  numerical  methods  for  solving  these  equations 
are  very  important  for  many  scientific  and  engineering  applications.  In  this 
paper  a  new  algorithm  is  presented  for  solving  finite  difference  equations  for 
the  linearized  Navier-Stokes  equations.  The  method  is  similar  to  successive- 
over-relaxation  which  is  a  widely  used  algorithm  for  solving  elliptic 
difference  equations.  Numerical  results  showing  the  behavior  of  the  method 
are  presented.  Other  results  appeared  in  an  earlier  report  which  discussed 
finite  difference  schemes  for  the  incompressible  Navier-Stokes  equations.  The 
method  is  efficient  and  easy  to  implement. 
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AM  ITERATIVE  METHOD  FOR  SOLVING  FINITE  DIFFERENCE 
APPROXIMATIONS  TO  THE  STOKES  EQUATIONS 

John  C.  Strikwerda 

1 .  Introduction 

In  this  paper  we  present  and  analyze  a  new  iterative  method  for  solving 
finite  difference  approximations  to  the  steady  Stokes  equations .  The  method 
is  a  variant  of  successive-over-relaxation  (S.O.R.)  and  has  similarities  to 
the  method  used  by  Chorin  (1968)  for  the  time-dependent  Navier-Stokes 
equations.  The  method  described  here  is  called  extended  successive-over- 
relaxation  (B. S.O.R.)  and  is  useful  for  solving  the  nonlinear  incompressible 
Navier-Stokes  equations  as  well. 

The  Stokes  equations  are 

-  v2  S+  %-  { 

(1.1)  in  a  c 

$  *  u  •  g 

and  we  take  as  boundary  conditions 

u  ■  b  on  9u  . 

♦ 

The  velocity  u  is  a  vector  of  dimension  k  and  the  pressure  p  is  a 
scalar.  The  system  (1.1)  requires  k  boundary  conditions  which  can  be  either 
of  Dirlchlet  type,  as  given  above,  or  some  other  type. 

A  commonly  used  method  for  solving  (1.1)  is  to  replace  the  second 
equation  of  (1.1),  the  divergence  equation,  by  an  elliptic  equation  for  the 
pressure.  The  resulting  finite  difference  approximation  can  be  solved  by 
iterative  methods  for  elliptic  equations,  (e.g.  Harlow  and  Welch  (1965), 

Roache  (1972)).  The  difficulty  with  this  approach  is  that  solutions  to  the 
derived  'system  need  not  be  solutions  of  the  original  system  (1.1)  (see 
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Strikwerda  (1983),  and  Greenspan  et  al.  (1964)).  Therefore,  we  consider  only 
finite  difference  approximations  to  the  system  (1.1)  in  the  form  given  there. 

The  finite  difference  approximation  of  (1.1)  results  in  matrix  equations 
of  the  form 


where  the  matrices  -Ah,  G^,  and  result  from  finite  difference 

approximations  of  the  vector  Laplacian,  gradient,  and  divergence  operators, 
respectively.  will  be  assumed  to  be  a  square  n  by  n  matrix,  G^  an 

n  by  art-1  matrix,  and  an  art-1  by  n  matrix.  He  will  denote  the 

n+mtl  by  n+nrt-1  matrix  in  (1.2)  by  Z^.  Systems  of  the  form  (1.2)  also 
arise  in  solving  the  time-dependent  Stokes  and  Navier-Stokes  equations. 

The  equations  (1.1)  will  not  have  a  solution  unless  the  integrability 
condition 

d.3)  Ja g  -  •  n 

is  satisfied.  Similarly,  the  matrix  Zh  in  (1.2)  will,  in  general,  be 
singular,  and  the  system  (1.2)  will  not  have  a  solution  unless  the  data  are 
constrained  to  be  orthogonal  to  the  left  null  vectors  of  the  matrix  Z^.  He 
will  assume  that  the  rank  deficiency  of  Z^  is  only  one,  corresponding  to  the 
one  integrability  condition  (1.3). 

Rather  than  constraining  the  data  to  be  orthogonal  to  the  left  null 
vector  of  Z^,  we  prefer  to  consider  both  p^  and  g^  as  defined  only  up  to 
arbitrary  additive  constants.  That  is,  they  are  elements  of  the  vector 
spaces  where  the  quotient  space  is  defined  by  v^  =  v2  if  v.,  -  v2 

has  all  components  equal.  Considered  this  way,  Z^  is  an  n+m  by  n+m  non¬ 
singular  matrix.  This  approach  to  solving  (1.1)  is  discussed  in  more  detail 
in  section  4  of  Strikwerda  (1983). 


The  E.S.O.R.  algorithm  discussed  in  this  paper  has  been  used  to  solve 
several  test  problems  involving  the  Stokes  equations  (8trikwerda  (1983))  and 
is  being  used  by  the  author  to  solve  for  solutions  of  both  the  steady  and  time 
dependent  Navier-Stokes  equations*  The  method  appears  to  be  quite 
efficient*  Dse  of  B.S*O.R.  as  a  pre-conditioner  for  a  conjugate  gradient 
algorithm  is  being  investigated* 

In  the  next  section  we  will  analyse  a  class  of  iterative  methods  for 
systems  of  the  form  (1*2)  and  in  section  3  we  will  discuss  how  the  methods 
behave  as  the  mesh  size  varies*  Several  numerical  examples  which  illustrate 
the  utility  of  the  method  are  given  in  section  4. 


2.  The  Extended  S.O.R.  Method. 

In  this  section  we  study  a  class  of  iterative  methods  to  solve  linear 
systems  of  the  form 


(2.1) 


where  X,  G,  and  D  are  matrices  of  dimension  n  x  n,  n  x  m,  and  m  x  n, 
respectively.  We  assume  that  the  matrix 


Z  - 


G 

0 


is  non-singular,  hence  n  is  greater  than  or  equal  to  m. 

Systems  of  the  form  (2.1)  often  arise  in  the  solution  of  constrained 
optimisation  problems.  Indeed,  the  solution  of  the  Stokes  equations  may  be 
regarded  as  the  minimum  of  a  quadratic  functional  under  the  constraint  that 
the  divergence  of  u  is  specified.  If  m,  the  number  of  constraints,  is 
much  smaller  than  n  one  can  often  eliminate  m  values  of  the  unknown  x 
using  the  second  row  of  equations  in  (2.1)  and  so  obtain  a  system  which  can  be 
solved  by  standard  methods,  see  e.g.  Dyn  and  Ferguson  (1982).  if  a  cartesian 
grid  is  used  on  a  rectangular  region  one  can  use  a  special  technique  developed 
by  Aait,  Hall,  and  Porsching  (1980)  to  reduce  the  system  to  an  n  by  n 
system,  we,  however,  will  consider  the  case  where  m  is  quite  large  and 
where  there  is  not  a  natural  or  convenient  way  to  reduce  the  system  to  one 
involving  only  n  equations  in  n  variables. 

The  iterative  methods  we  will  discuss  are  extensions  of  successive-over- 
relaxation  (S.O.R.)  as  applied  to  the  matrix  A.  We  assume  that  X  has  been 
transformed  so  that 


X  -  I  -  I.  -  U 


where  L  and  U  are  strictly  lower  and  upper  triangular  matrices, 
respectively*  The  S.O.R*  iterative  procedure  applied  to  the  system 


Ax 


Vfl  V  ,  V  .  Vfl  V  . 

x  -  x  -  u(x  -  Lx  -  Ox  -  a) 


is  given  by 
(2*2) 

and  we  assume  that  this  converges  for  u  satisfying  0  <  w  <  for  some 
positive  value  of  For  the  basic  theory  of  S.O.R.  the  reader  is  referred 

to  Young  (1971).  Assuming  that  (2.2)  converges  is  equivalent  to  assuming  the 
following  condition. 

Condition  2.1 


There  is  a  positive  constant  aig  such  that  for  0  <  w  <  the  roots  of 

(2.3)  det(— — 1  I  -  XL  -  U)  -  0 

u 

satisfy  | X|  <  1. 


For  the  full  system  (2.1)  we  consider  the  extended  successive-over¬ 
relaxation  iterative  procedure 


vfl  v  ,  v  vfl  _  v  .  _  v  .  _  v+1 

x  ■  x  -  w(x  -  Lx  -  Ox  f  Gfty  +  G^y  -  a) 


(2.4) 


where 


vfl  v  v  .  _  vfl  .. 

T  -  Y  -  Y(Dq*  +  DjX  -  b) 


and 


G0  +  Gt 


Oq  f  D1  -  D  . 

The  iterative  parameters  must  be  determined  so  that  (2.4)  is  a  convergent 
algorithm.  The  purpose  of  the  analysis  in  this  section  is  to  find  conditions 
under  which  (2.4)  will  converge. 


Chorin  (1968)  used  a  scheme  similar  to  (2.4)  to  solve  for  the  velocity 
and  pressure  at  each  new  time  level  for  the  time- dependent  Navier-Stokes 
equations.  In  Chorin' s  method  the  matrix  A  is  essentially  the  identity 
matrix  and  he  set  u>  to  be  1.0  and  and  Dg  were  zero, 

lie  rewrite  (2.4)  in  matrix  form  as 


(2.5) 


where 


V  +  c 


1-UI 

—  1+0  ^0 


The  method  (2.5)  will  converge  if  and  only  if  all  the  eigenvalues  of  X^  XQ 
have  absolute  value  less  than  one.  The  first  result  on  the  eigenvalues  of 
X^1X0  is  this  lemma. 


For  y  ■  0  there  are  two  classes  of  eigenvalues  of  X^  X^.  There  are 
n  eigenvalues  which  are  roots  of  (2.3)  and  m  simple  eigenvalues  all  equal 
to  unity. 

Proof 


Let  X  be  an  eigenvalue  of  x^'x^,  then 

0  -  detUl  -  x7 1  X„ )  -  det(  Xx„  -  X.  )det  xT1 


■o  X1  is  non-singular  for  small  values  of  y*  We  have 


(2.6) 


*1  -xo 


X+arl 


I  -  XL  -  0  X0t  +  Gq 


Y(XD1  ♦  Dq) 


(X-1)I 


so  at  Y  ■  0  the  eigenvalues  of  X^  XQ  are  either  roots  of  (2*3)  or  are 
equal  to  unity.  The  eigenvalues  equal  to  unity  are  easily  seen  to  be  simple 
because  A  is  non-singular.  In  fact/  for  any  y  e  rf1  the  vector 


is  an  eigenvector  of  X^  XQ  at  Y  “  0.  This  proves  Tamms  2.1. 

The  n  eigenvalues  of  X~^Xg  which  are  equal  to  the  roots  of  (2.3)  at 
Y  “  0  will  be  called  the  8.0. R.  eigenvalues  of  x^Xg. 

We  will  now  study  the  perturbation  expansion  of  the  eigenvalues  of 


Xj  Xg  about  Y  *  0< 
Theorem  2.1 


y 


a 


Those  eigenvalues  X^(y)  of  X^  XQ  such  that  X^(0)  -  1  satisfy 
(2.7)  X^y)  -  1  -  ly  ♦  0(Y> 

where  is  an  eigenvalue  of  -DA-1G. 

The  proof  of  Theorem  2.1  depends  on  the  following  two  lessaas. 


Let  T(y)  be  an  analytic  matrix-valued  function  defined  in  a  neighbor¬ 
hood  of  y  ■  0.  If  Xq  is  a  simple  eigenvalue  of  T(0)  then  the  eigen¬ 
values  X^(y)  of  T(y)  for  which  X^(0)  *  Xg  satisfy 

"  \>  +  YMj  *  °(y)  ' 

where  is  at.  ,.«nv  am  at  T'(0). 


>vv;v. 


.v.-. .%  % v. 


Lemma  2.3 


If  a,  b,  c,  and  d  are  matrices  of  dimension  n  x  n,  n  x  m,  m  x  n. 


and  m  x  m,  respectively,  then 


det(a  . }  =  det(a  -  bd  c)  det  d,  if  det  d  ?  0 


det(d  -  ca  b)  det  a,  if  det  a  j*  0 


Proof  of  Lemma  2.2. 

The  proof  easily  follows  from  results  of  Kato  (1966)  but  we  give  it  here 


for  completeness .  For  y  near  zero  we  can  find  a  non-singular  analytic 
matrix  valued  function  P(y)  such  that 


p(y)  t(y)p(y)  *  t(y) 


tq(y)  o 


has  a  block  form 


\  0  T^y)/ 

where  T(0)  -  XQI.  We  now  consider  only  Tq(y).  The  eigenvalues  of  Tq(y) 


have  expansions  as  Puiseux  series 


YTl '  xo  *  J,  v*/p 

for  some  positive  integer  p.  If  p  ■  1,  the  result  follows.  Assume  that 
p  >  1,  we  will  show  that  for  1  <  4  <  p,  \  f  is  zero. 


V*’  ‘  U3*',‘/P  '  UJ0  *  0 


be  the  eigenvector  corresponding  to  X^(Y>*  (u^ ( y)  does  not  have  a  pole  at 

Y  ■  0  since  TQ(0)  is  diagonal.)  Since  Tq(y)  “  XQI  +  YTQ1  +  Y2TQ2  +  ••• 


T„<Y)u,(y)  -  uj(y)Xj(y)  , 


Since  the  iterative  procedure  (2.4)  will  be  convergent  only  if  the 
eigenvalues  of  x^xq  are  all  less  than  one  in  magnitude,  we  see  from  Theorem 
2.1  that  for  (2.4)  to  converge  for  positive  values  of  y  we  must  assume  that 
the  following  condition  holds. 

Condition  2.2. 

All  the  eigenvalues  of  -DA-^G  have  positive  real  part. 

Mote  that  if  all  eigenvalues  of  -DA- ^G  have  negative  real  part,  then 
one  can  either  multiply  the  last  m  equation  of  (2.1)  by  negative  one  (i.e. 
replace  D  by  -to)  or,  equivalently  take  Y  to  be  negative.  If,  however, 
some  of  the  eigenvalues  of  -DA-1G  have  positive  real  part  and  others  have 
negative  real  part  then  the  method  will  not  converge. 

He  now  state  the  main  result  of  this  section. 

Theorem  2.3 

Conditions  2.1  and  2.2  are  sufficient  for  the  algorithm  (2.4)  to  converge 
for  Y  and  w  satisfying  0  <  Y  <  Y„  and  0  <  u  <  uQ  for  some  positive 
values  of  Yq  and  Furthermore,  if  A  is  non-singular  then  necessary 

conditions  for  (2.4)  to  converge  for  such  f  and  u  are  that  the  roots  of 
(2.3)  satisfy  |1|  <1  and  that  the  eigenvalues  of  -DA-1G  have  non-negative 
real  part. 

Proof 

Consider  the  two  groups  of  eigenvalues  described  in  T<enma  2.1.  By 
continuity  of  the  eigenvalues  as  functions  of  the  matrix  elements  we  have  that 
the  8.O.R.  eigenvalues  satisfy  | A|  <  1  for  y  in  some  range  0  <  Y  <  Yq 
for  0  <  u  <  <i)q.  Then  by  Theorem  2.1  and  Condition  2.2  the  remaining  eigen¬ 
values  of  also  satisfy  1 1|  <  1  for  0  <  Y  <  YQ  for  some  yq*  This 


proves  the  sufficiency  condition 


then 


If  the  algorithm  (2.4)  converges  for  0  <  Y  <  Yq  and  0  <  at  <  u>Q, 
for  Y  *  0  the  eigenvalues  of  x^xq  must  satisfy  |  X |  <  1.  Since  A  is 
non-singular  the  S.O.R.  eigenvalues  of  x^Xq  are  not  equal  to  1  for 
Y  ■  0.  Thus  the  non-S.O.R.  eigenvalues  are  simple  and  satisfy  (2.7).  The 
condition  that  |X|  <  1  for  small  positive  y  implies  that  Re  tk  >0. 

We  conclude  this  section  by  obtaining  expressions  for  the  perturbation  of 
the  S.O.R.  eigenvalues  of  the  iteration  matrix  for  the  case  where  A  is 
diagonal! sable  and  has  property  A.  Under  these  conditions  the  iteration 
matrix  for  S.O.R.  applied  to  A  has  principle  vectors  of  grade  at  most  two, 
(Young  (1971)  p.  233-238). 

Now  let  be  a  simple  S.O.R.  eigenvalue  for  x^x0  at  y  m  0  and 


be  the  perturbation  expansion  of  a  right  eigenvector  of  X^Xg  with  eigen¬ 
value 


(2.8)  -  A^0  +  Y^jj  +  0(Y)  • 


be  a  left  eigenvector  at  y  m  0  such  that 


(vj0'uj0)  ?  0  * 

Substituting  the  above  expansions  for  the  eigenvector  and  eigenvalue  in  the 
equation 


we  obtain 


(Wj,XqWj)  ■  Ajtw^XjWj) 


x  m  _  <vjo'(ywpii) 
31 


11- 


3.  The  Finite  Difference  Stokes  Equations. 

In  this  section  we  consider  the  application  of  the  iterative  method  (2.4) 
to  finite  difference  approximations  of  the  Stokes  equations.  We  first 
consider  Conditions  2.1  and  2.2  to  see  if  they  are  satisfied. 

Since  the  matrix  A  arises  from  a  discretization  of  the  vector 
Laplacian,  Condition  2.1  is  very  reasonable.  If  the  finite  difference  grid  is 
rectangular  with  uniform  spacing  and  one  uses  the  standard  five-point 
discretization  for  the  Laplacian,  then  Condition  2.1  is  satisfied.  Young 
(1971).  In  addition,  A  will  be  symmetric  and  have  Property  A,  (Young 
(1971)). 

Condition  2.2  will  also  be  satisfied  for  appropriate  difference 

schemes.  The  operator  Qjj  represented  by  is  a  f*nite  difference 

2 

approximation  to  the  operator  Q0  defined  on  L  ( fi)/R  as  follows. 

Q0P  -  q  if 
q  =  $  •  u 

where  V2  u  =  $p  in  ft 

with  u  =  0  on  3ft  . 

Crozier  (1974)  has  proved  the  following: 

Theorem  3.1 

If  SI  is  a  connected,  bounded  domain  in  V?  with  smooth  boundary  then 

2 

the  operator  Qg  is  a  bounded,  positive  definite  operator  on  L  (Q)/R. 

Therefore,  if  is  a  consistent  approximation  to  Qg  one  can  expect 

that  the  next  condition  holds. 

Condition  3.1. 

There  are  positive  constants  c.,  and  c2  such  that  for  0  <  h  <  hQ, 

(3.1)  c1  <  Re  rij  <  InJ  <  c2 

where  the  n.  are  as  in  Theorem  2.1. 


It  is  important  to  note  that  Condition  3.1  is  not  satisfied  for  all 


finite  difference  schemes.  In  particular,  if  one  uses  standard  central 
difference  to  approximate  both  the  gradient  of  the  pressure  and  the  divergence 
of  the  velocity,  then  numerical  tests  indicate  that  Condition  3.1  is  not 
satisfied.  In  section  4,  we  will  discuss  difference  schemes  which  satisfy 
Condition  3.1.  Condition  3.1  is  related  to  the  regularity  of  the  difference 
scheme  (Bube  and  Strikwerda,  1983).  Regular  difference  schemes  are  those 
whose  solutions  satisfy  regularity  estimates  analogous  to  those  satisfied  by 
solutions  of  the  differential  equation. 

We  now  consider  the  convergence  behavior  of  the  K. S.O.R.  method.  Suppose 
then  that  one  has  a  finite  difference  approximation  to  the  Stokes  equations 

(1.2)  for  which  the  method  (2.6)  will  converge  for  sqm  positive  values  of  y 
and  u.  One  would  like  to  know  how  to  choose  values  of  m  and  y  so  as  to 
obtain  a  good  rate  of  convergence  for  the  method.  We  are  unable  to  give 
rigorous  estimates  of  the  convergence  rate,  but  we  will  now  show  that  the 
following  conjecture  is  quite  plausible. 

Conjecture  3.1 

If  the  matrix  K  satisfies  property  h  and  Condition  3.1  is  satisfied 
then  there  are  positive  constants  cQ  and  c1  such  that  for  u  -  2/(1+cQh) 
and  y  »  c^h  then 

(3.2)  ptX^Xjj)  -  1  -  Kh  +  0(h) 
for  some  positive  constant  K. 

Since  the  iteration  matrix  for  S.O.R.  applied  to  the  discrete  five  point 
Laplacian  satisfies  a  relation  like  (3.2),  Conjecture  3.1,  if  true,  shows  that 
B. S.O.R.  for  the  Stokes  equations  is  roughly  as  efficient  as  S.O.R.  for  the 
five  point  Laplacian. 


He  argue  for  the  Conjecture  3.1  as  follows.  If  A  satisfies  Property  A 
•  * 

then  for  w  >  <■>  ,  where  w  is  the  optimal  parameter  for  (2.2),  and  small 
positive  values  of  y,  the  S.O.R.  eigenvalues  satisfy  (2.8)  or  (2.10)  and 
have  modulus  <*>-1  +  0(  y) .  Consider  first  those  A^Q  which  are  near  u  -  1, 


that  is,  XjQ  has  the  form 


Aj0  »  (<v-1)e 


where  8^  -  0(1)  and  u  ■  2  ^  0(h)  as  h  tends  to  zero.  If  the  finite 
difference  approximations  for  the  divergence  and  gradient  are  consistent  then 
(2.9)  can  be  approximated  by 


X  ** 

jl 


'v*1 

(Uj,  A  U^) 


since  the  adjoint  of  the  gradient  is  the  negative  of  the  divergence  and  where 
A  is  the  finite  difference  negative  Laplacian,  i.e.  without  the  normal¬ 
ization  which  makes  the  diagonal  elements  unity.  The  discrete  eigenvector 
Uj  may  be  regarded  as  a  representation  of  smooth  vector  function  u  and  so 
the  above  ratio  is  approximated  as 

x  *  l*1*  «l2 

| grad  $|2 

and  so  A^  is  0(1)  as  h  tends  to  zero. 

On  the  other  extreme  where  is  close  to  -(u-1)  the  discrete 

eigenvector  Uj  is  a  very  oscillatory  function.  Then  we  have 

(vyvuj0  -°(h’1) 

(VAjoGi)  • 0lh) 


AUjQ  -  0(1) 

and  thus  A^  is  again  0(1)  as  h  tends  to  zero. 
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For  other  values  of  X^  on  the  circle  with  radius  <u  -  1,  an  argument 

similar  to  those  above  shows  that  X^  will  be  bounded  as  h  tends  to 

zero.  For  X .  *,  as  in  (2.10)  and  (2.11),  the  conclusion  is  that  X.  *.  is 
j#/2  3  *'2 

V, 

proportional  to  h  ^  as  h  tends  to  zero. 

Therefore  if  y  is  taken  proportional  to  h  in  (2.8)  and  (2.10)  and  we 
assume  that  the  terms  which  are  0{y)  in  (2.8)  and  0( y)  in  (2.10)  become 
0(h)  and  0(h),  respectively,  we  then  obtain  (3.2).  This  last  assumption  is 
the  one  for  which  we  have  no  theoretical  justification.  It  is  not 


unreasonable,  however,  and  the  numerical  experiments  confirm  that  Conjecture 


3.1  is  quite  plausible 


4.  Numerical  Examples 


in  this  section  we  present  some  numerical  results  of  using  the  E.S.O.R. 
algorithm  on  a  test  problem.  We  consider  the  Stokes  equations 


V^u  -  &  m  0 
3x 


(4.1) 


V*V  -  -  0 


■|~  +  -rp  ■  g(x,y)  ”  cos  s  x  cos  i  y 


on  0  <  x,y  <  1  with  u  and  v  specified  on  the  boundary.  The  exact 
solution  is  given  by 


-1 


(2s)  sin  s  x  cos  ir  y 


.-1 


(2s)  cos  s  x  sin  s  y 


cos  s  x  cos  s  y 


The  discretization  used  a  uniform  grid  with  the  same  number  of  grid 
points  in  each  direction.  The  second-order  accurate  five-point  Laplacian  was 
used  to  approximate  the  Laplacian.  As  mentioned  in  section  3  the  choice  of 
discretization  for  the  gradient  and  divergence  terms  is  crucial  to  satisfy 
Condition  3.1.  We  employed  here  the  regularized  centered  differences 


(Strikwerda  (1983))  given  by 

JE 


dx  *  ~x0 


.  1.2.  .2 

o  _p  -  —  h  6  a  ,p 

vOr  ft  v— 


39  “  r  h*  4  p 

3y  yO*  6  y-  y+ 


3u  .  1  .2  .  .2 

— —  «  o  u  -  —  h  6  6  u 
3x  xO  6  x+  x- 


3v  .  1  w2  .  .2 

-r-*o.v--h  6  6  v 

3y  yO  6  y+  y- 


where  h  is  the  grid  spacing  and  5^ ,  6^,  and  6^_  are  the  centered, 


x+' 


forward,  and  backward  difference  operators  in  the  x-direction.  The  operators 


4  5  . ,  and  6  are  defined  similarly  for  the  y-direction.  To  determine 

yO  y+  y- 

the  pressure  at  boundary  points  a  cubic  interpolation  was  used,  e.g. 

Plj  =  3(P2j  “  P3j>  "  P4j  * 

In  these  tests  the  difference  operators  DQ  and  were  zero,  i.e. 

*  D,  Gq  *  G.  This  is  perhaps  the  easiest  scheme  to  implement  of  those 
considered  here.  Mote  that  for  the  theory  developed  in  section  2  it  is 
necessary  that  the  difference  operators  GQ  and  G^,  together  with  G, 
annihilate  constants  so  that  they  are  defined  on  H®*  V*.  This  scheme  has 
been  shorn  to  be  second-order  accurate,  Strikwerda  (1983).  Here  we  give 
results  only  on  the  efficiency  of  the  solution  algorithm,  ESOR. 

To  support  the  conjecture  3.1,  several  runs  were  made  where  <■>  and  y 
were  given  by 


(4.1) 

for  several  values  of 
the  quantities 


(4.2) 


u  -  2/(1  +  cQh) 

Y  "  eft 

Cq,  Cj  and  h.  The  iterative  method  was  stopped  when 

.  n+1  n,  .  n+1,2  , 

lu  -  u  I  /  (lu  I  +1)* 

.  n+1  n,  .  ,,  n+1 .2  .  *A 
Iv  -  v  I  /  (Iv  1  +  1)  *■ 

.n+1  n.  .  n+1 .2  . 

Ip  -  p  I  /  { Ip  I  +  1)  * 


were  all  less  than  10~^.  The  norms  for  u  and  v  in  (4.2)  were  the 

2  2 
discrete  L  norms,  and  the  norm  for  p  was  the  L  norm  in  the  quotient 

space  rf"+  /«.  The  computation  of  the  norm  in  the  quotient  space  will  be 

discussed  later.  If  Conjecture  3.1  were  valid  then  the  product  of  I,  the 

number  of  iterations  required  for  convergence,  and  h  would  tend  to  a  limit 

as  h  tends  to  zero.  To  see  this  we  observe  that  if  the  spectral  radius  is 

1  -  Kh  then  the  number  of  iterations  required  for  the  relative  change  in 


successive  iterates  to  he  less  than  e  is  determined  by 


(1  -  Kh)  »  e 


This  implies 

(4.3)  hi  »  -(log  e)/K  . 

The  results  of  these  runs  are  shown  in  Table  1 . 

The  results  in  Table  1  for  cQ  =  5.0  and  c1  -  5.0  give  excellent 
agreement  with  Conjecture  3.1.  The  variation  in  the  values  of  h«I  for  other 
values  of  Cq  and  c^  can  be  explained  by  the  presence  of  the  o(h)  term  in 
(3.2)  and  because  the  use  of  the  norms  is  only  an  imperfect  indicator  of  the 
spectral  radius. 

We  now  discuss  the  computation  of  the  norms  for  the  quotient  spaces 
If  X  is  a  vector  in  then  the  L2  norm  of  its  image  in 

if*1/*  is 

“+1  _  2  V2 

IX!  -  (  l  (X  -  xr) 
k«1 

where  X  is  the  mean  of  X,  i.e. 


.  m+1 

x  “  Sm  i  Xk  * 

k-1 

An  efficient  and  accurate  way  to  compute  this  norm  is  the  algorithm  of  West 
(1979).  This  can  be  described  as  follows 
Initialise  k  =  1 


repeat  for 


*1,-0 


k  —  2, , m+ 1 
R  -  xk  - 


U 


rA 


"k-i 


+  u 


Nk  *  Nk-1  +  R  *  u  *  (k-1) 


finally 


4X1  -  (N 


>  . 


This  algorithm  is  stable  as  shown  by  Chan  and  Lewis  (1979),  and  is  very 
convenient.  This  is  used  to  compute  the  norm  of  the  pressure  and  the  residual 
of  the  last  equation  in  the  Stokes  equations. 

The  mean  of  the  residual  of  the  divergence  equation  is  the  quantity  6^ 
discussed  in  Strikwerda  (1983),  and  for  each  of  the  cases  reported  here  it  was 
on  the  order  of  the  truncation  error. 

The  E.S.O.R.  method  described  here  was  used  to  compute  the  solutions 
discussed  in  Strikwerda  (1983)  where  accurate  finite  difference  schemes  for 
the  Stokes  equations  are  described.  Xt  is  also  being  used  in  the  computation 
of  Taylor  vortex  solutions  to  the  steady  Navier-Stokes  equations  and  in 
computations  of  solutions  of  the  time-dependent  Navier-Stokes  equations.  The 
results  of  this  research  will  be  reported  when  it  is  completed,  it  has  been 
found  to  be  a  reliable  algorithm,  the  main  difficulty  being  the  choice  of 
<■>  and  y.  Conjecture  3.1  provides  a  means  of  estimating  good  values  of 
o>  and  y.  By  finding  good  values  of  w  and  y  for,  say,  h  *  1/10,  one 
can  then  use  Conjecture  3.1  to  obtain  good  estimates  of  u  and  y  for 


smaller  values  of  h 


5.  Conclusion 


The  E.S.O.R.  method  has  been  rigorously  analyzed  for  matrices  of  the  form 
(2.1)  with  the  main  results  stated  in  Theorem  2.3.  For  the  particular  case  of 
difference  approximations  to  the  Stokes  equations  we  have  argued  in  section  3 
that  the  assusiptions  required  by  Theorem  2.3  are  reasonable  for  many  finite 
difference  schesws.  The  results  of  section  4  have  confirmed  that  the  E.S.O.R. 
method  is  indeed  an  efficient  algorithm  for  the  solution  of  the  finite 
difference  Stokes  equations. 
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